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ABSTRACT: We perform an extensive study of autocorrelation of several observables in lattice 
QCD with two degenerate flavours of naive Wilson fermions and unimproved Wilson gauge action 
using DD-HMC algorithm. We show that (1) at a given lattice spacing, autocorrelation of topo- 
logical susceptibility decreases with decreasing quark mass and autocorrelations of plaquette and 
(y-) Wilson loop do not increase with decreasing quark mass, (2) autocorrelation of topological sus- 

ceptibility substantially increases with decreasing lattice spacing but autocorrelation of topological 
charge density correlator shows only mild increase and (3) increasing the size and the smearing 
level increase the autocorrelation of Wilson loop. 
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1 Introduction 

The most popular algorithm to simulate lattice QCD with Dynamical fermions is the Hybrid Monte 
Carlo (HMC) [1] and one of its improved variations, namely, Domain Decomposed Hybrid Monte 
Carlo (DD-HMC) [2] aims to achieve significant acceleration of the numerical simulation. Dynam- 
ical Wilson fermion simulations at smaller quark masses, smaller lattice spacings and larger lattice 
volumes on currently available computers have become feasible with recent developments such as 
DD-HMC algorithm. However, approach to the continuum and chiral limits may still be hampered 
by the phenomenon of critical slowing down. One of the manifestation of critical slowing down 
is the increase in autocorrelation times associated with the measurements of various observables. 
Thus measurements of autocorrelation times help us to evaluate the performance of an algorithm 
in terms of critical slowing down. In addition, an accurate determination of the uncertainty asso- 
ciated with the measurement of an observable requires a realistic estimation of the autocorrelation 
of the observable which in turn depends on the various parameters associated with the particular 
algorithm used. 

An extensive study of autocorrelation mainly in pure St/ (3) gauge theory using DD-HMC 
algorithm has been carried out by ALPHA collaboration [3]. They have shown that the autocor- 
relation of squared topological charge increases dramatically with decreasing lattice spacing while 
Wilson loops decouple from the modes which slow down the topological charge as lattice spacing 
decreases. In the simulations with dynamical fermions, the study becomes more difficult, because 
the autocorrelation may now depend on number of quark flavours (tif), the quark masses and the 
fermion action used [4]. In fact ALPHA collaboration [3] has shown, in the case of tif = 2 QCD 
with Clover action for a given value of quark mass and lattice volume, that squared topological 
charge decorrelates faster compared with pure gauge at approximately same lattice spacing. These 
dependencies and the one on the lattice spacing remain to be studied in detail. In this work we study 
the autocorrelations of a variety of observables measured with DD-HMC algorithm in the case of 
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naive Wilson fermions [5,6]. Note that the measurement of autocorrelation is notoriously difficult, 
since accurate determination of it may require considerably larger accumulated statistics (total 
molecular dyanamics time). In this work we mainly focus on various trends of autocorrelations 
that we can observe clearly rather than the precise measurement of the integrated autocorrelation 
time. 



2 Autocorrelation 



Following Refs. [3] and [8] let us assume that & = | &(x) \ be a real-valued function defined 
on the state space 5 that is square integrable with respect to %, where % is the stationary Markov 
chain probability distribution with probability transition matrix P. Now consider that the Markov 
chain is in equilibrium. Then the unnormalized autocorrelation function, 



c G {t) = {&{ S ms+t))-^ 

1 0(y) (2.1) 



AM' 



ftxPxy Tt x 1t y 



where = {&{t)) = Ya K x &(x). Now if the algorithm satisfies detailed balance, i.e. K x P xy = 7i y P yx 
for all x,y £ 5 then it is convenient to introduce the symmetric matrix 

T x , y = 7ihx y n y ~ i[ (2.2) 

which has real eigenvalues X n , n > with Ao = 1 and | Xn |< 1 for n > 1, assuming an ergodic 

algorithm. We order the eigenvalues as | X n \<\ A„_i |. There is a complete set of eigenf unctions 

i 

%n{x) with %o{x) = 7i x . By using spectral representation of T, Eq. (2.1) can be reduced to 

C^) = £(A„)' |T] n (^| 2 (2.3) 

n>\ 

1 

where r\ n {6)= £ x & \x)x n (x)% x . Since | Xn \< 1 for n > 1 

<^(0 = IX' /T "l^(^)| 2 (2-4) 

n>l 

where T„ = — j4r-, assuming A„'s are positive. 

For any particular observable G, autocorrelation among the generated configurations are gen- 
erally determined by the integrated autocorrelation time T^ t for that observable. For this purpose, 
at first, one needs to calculate the unnormalized autocorrelation function of the observable & mea- 
sured on a sequence of ,/V equilibrated configurations as 

C e (0 = - J- £ {ff r - G) {@ r+t - W) (2.5) 

ly 1 r=\ 

where & = A YH-=\ &r is the ensemble average. Following the windowing method as recommended 
by Ref. [8], the integrated autocorrelation time is defined as 

i w 

< = \ + ^{t) (2.6) 
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where T (?) = C (?) jC (0) is the normalized autocorrelation function and W is the summation 
window. To calculate the errors, we follow the standard techniques available in the literature [4, 8- 
11]. The variance of T & (?) is given by 

1 °° 

((8F e (t)) 2 } « - £ [T e (k + t)+ T e (k-t)- 2T e (?) T e (k)} 2 (2.7) 
and the variance of %f nV 

/fst r ff\2\ „ 2(2W + 1) e 2 

\(<> T mf)>~ K^int) ■ ( 2 - 8 ) 

Different strategies have been suggested in the literature [4, 8, 11] for choosing W. We choose W 
where error of (t) becomes equal to r^(?) [4]. The above expressions are used to calculate 
the errors unless otherwise stated. In case the total accumulated statistics is extremely large an 
alternative procedure may be to use binning, with binsizes much larger than ii nt for calculating the 
error [7]. 



3 Observables 

Let us denote plaquette and Wilson loop of size RxT with smear level s by P s and W S (R,T) 
respectively. Topological susceptibility with smear level s is denoted by Q 2 (the normalization fac- 
tor, inverse of lattice volume, is ignored). We have measured the autocorrelations for the plaquette, 
Wilson loop, nucleon propagator, pion propagator, topological susceptibility and topological charge 
density correlator (C(r) = (q(x)q(0)) where q(x) is topological charge density and r =| x |) for the 
saved configurations except for the unsmeared plaquette where we have measured for all the con- 
figurations, at two values of gauge coupling (j6 = 5.6 and 5.8) and several values of the hopping 
parameter K. 

For pion and nucleon we consider the following zero spatial momentum correlation functions 

d{t) = (01^(0^(0)10) and C 2 (t) = <0 | ^?(/)<7 2 (0) | 0) (3.1) 

where t refers to Euclidean time. For the nucleon O = N^N with ,/V = {q T d Cy^q u )q u . For the 
pion G^G = P^P or A 1 " A and (^i) 1 ' ff 2 =A^P or P f A with P = qflsqj and A 4 = qj^qj denote 
the pseudoscalar density and fourth component of the axial vector current (i and j stand for flavor 
indices for the u and d quarks and for the charged pion i ^ j). For both pion and nucleon we use wall 
source and point sink. We measure the autocorrelation of the zero spatial momentum correlation 
functions at an appropriate time slice corresponding to the plateau region of the effective mass. For 
lattice volume 24 3 x 48 and 32 3 x 64 we use 12 th and 15 th time slices respectively. For topological 
charge density, we use the lattice approximation developed for SU (2) by DeGrand, Hasenfratz and 
Kovacs [12], modified for St/ (3) by Hasenfratz and Meter [13] and implemented in the MILC 
code [14]. To suppress the ultraviolet lattice artifacts, smearing of link fields is employed. Unless 
otherwise stated 20 HYP smearing steps with optimized smearing coefficients a = 0.75, a 2 = 0.6 
and CC3 = 0.3 [15] are used for the gauge observables. For observables with hadronic operators 
no gauge field smearing has been used. Our data for topological charge, susceptibility and charge 
density correlator are presented in [16-18] 
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Table 1. Lattice parameters and simulation statistics. Here block, N2, N tr j, T refers to DD-HMC block, step 
number for the force F2, number of DD-HMC trajectories and the Molecular Dynamics trajectory length 
respectively. 

4 Auto-correlation Measurements 

We have generated ensembles of gauge configurations by means of DD-HMC [2] algorithm using 
unimproved Wilson fermion and gauge actions [5, 6] with tif = 2 mass degenerate quark flavors. 
At j6 = 5.6 the lattice volumes are 24 3 x 48 and 32 3 x 64 and the renormalized physical quark mass 
(calculated using axial Ward identity) ranges between 25 to 125 MeV (MS scheme at 2 GeV). 
At j6 = 5.8 the lattice volume is 32 3 x 64 and the renormalized physical quark mass ranges from 
15 to 75 MeV. To determine the lattice spacing we plot the ratio of lattice pion mass to lattice 
nucleon mass versus lattice pion mass. Extrapolation of the ratio to the physical point gives the 
lattice spacing a. The lattice spacings at jS =5.6 and 5.8 are 0.069 and 0.053 fm respectively. The 
Sommer method of determining the scale agree with the quoted value of lattice spacings at j3 = 5.6 
and P = 5.8 for the value of Sommer parameter ro = 0.44 fm. 

The number of thermalized configurations ranges from 7000 to 14000. The lattice parameters 
and simulation statistics are given in Table 1. For all ensembles of configurations the average 
Metropolis acceptance rates range between 75 — 98%. 

5 Results 

In Fig. 1 we show the autocorrelation function for the unsmeared plaquette for the total accu- 
mulated statistics 500, 1000 and 5620 respectively for j8=5.6, K = 0.15825 and lattice volume 
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Figure 1. Autocorrelation functions for the unsmeared plaquette for the total accumulated statistics 500, 
1000 and 5620 respectively at j3 = 5.6 for the ensemble B^j,- 

24 3 x 48. We notice that for smaller statistics, the autocorrelation function touches zero earlier 
leading to the underestimation of Xt nt . Also the positivity of the autocorrelation function is violated 
in contrast to theoretical expectations but the situation improves as statistics increases. 

Since it is exorbitant to measure smeared Wilson loops, propagators and smeared topological 
charge on each and every trajectory, we have measured these observables for the configurations 
saved with specific gaps. However unsmeared plaquette (Po) is measured on each trajectory. It 
is mandatory, however, to check that the measured autocorrelation scales appropriately with the 
gaps, so as to ensure the correct determination of the autocorrelation. We have carried out such 
checks and a typical result is presented in Figs. 2. In Fig. 2 we present the normalized autocorrela- 
tion functions (left) and integrated autocorrelation times (right) for unsmeared plaquette with three 
different gaps (4,8,32) between measurements for the ensemble C5. The data clearly exhibit the 
scaling properties with the gaps. 

In Figs. 3 and 4 we show noramlized autocorrelation functions and integrated autocorrelation 
times for Q\ Q at different fc's for j8 = 5.8 and j8 = 5.6 respectively. Windows are chosen as in- 
dicated by the vertical lines. Figs. 3 and 4 show that at both the lattice spacings (/3 = 5.6,5.8) 
autocorrelations of Q\ Q decrease with decreasing quark mass even though for the smaller quark 
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Figure 2. Normalized autocorrelation functions (left) and integrated autocorrelation times (right) for un- 
smeared plaquette with three different gaps (4, 8, 32) between measurements for the ensemble C5 at j3 = 5.6. 




Figure 3. Normalized autocorrelation functions (left) and integrated autocorrelation times (right) for at 
/3 = 5.8 for the ensembles D\ and D$. 

mass at j8 = 5.8 (D5) molecular dynamics trajectory length (t) is smaller. Note that the trend is 
more visible at smaller lattice spacing (/3 = 5.8). A possible explanation 1 for this suppression 
of autocorrelation with decreasing quark mass is that the algorithm needs to span between lesser 
number of topological sectors at smaller quark mass since the width of the Gaussian distribution of 
topological charge decreases with decreasing quark mass. 

In Fig. 5 we show normalized autocorrelation functions for C{r = 12) and Q\ Q at j3 = 5.6 for 
the ensemble B^y. Fig. 5 shows that at jS = 5.6 the autocorrelations for Q\ Q and C(r = 12) are very 
close. In Fig. 6 we show normalized autocorrelation functions and integrated autocorrelation times 
for C(r = 12) (left) and (right) at /3 = 5.8 for the ensemble D\ where pion mass is comparable 
with the pion mass for the ensemble B^- Fig. 6 shows that at j8 = 5.8 the autocorrelation for Q\ 
is larger than the autocorrelation for C(r =12). Figs. 5 and 6 show that autocorrelation for Q\ Q 
increases quite significantly with decreasing lattice spacing at comparable pion mass whereas the 
autocorrelation of topological charge density correlator (C(r)) increases slightly with decreasing 

1 Stefan Schaefer (private communication). 
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W= 180 





Figure 4. Normalized autocorrelation functions and integrated autocorrelation times for (?2o at P =5.6 for 
the ensembles By, (left), (right) and B^ (bottom). 



lattice spacing. Taking into account the effect of active link ratio (see for example section 3.1 in 
Ref. [3]) (R = 0.363 for j8 = 5.6 and R = 0.422 for j8 = 5.8) strengthens this conclusion. 

In Fig. 7 we present normalized autocorrelation functions (left) and integrated autocorrelation 
times (right) for Pq for j8 = 5.6 at lattice volume 24 3 x 48. We find that Ti nt for Pq is not increasing 
with decreasing quark mass. In Fig. 8 we present normalized autocorrelation functions (left) and 
integrated autocorrelation times (right) for W^o(4,4) for j8 = 5.6 at lattice volume 24 3 x 48. The 
figure shows T,„ ? for W2o(4,4) is also not increasing with decreasing quark mass. 

For the measurement of static potential V(r) one needs to measure Wilson loops of various 
sizes. In the measurement of a Wilson loop, to suppress unwanted fluctuations smearing is needed. 
Therefore it is interesting to study how autocorrelation of Wilson loops changes with sizes of the 
Wilson loops and smearing levels. In Fig. 9 we present normalized autocorrelation functions and 
integrated autocorrelation times for W20 with different sizes for the ensemble Dy. In Fig. 10 we 
show normalized autocorrelation functions and integrated autocorrelation times for W(3,3) with 
different levels of HYP smearing for the ensemble C2. We observe that the integrated autocor- 
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Figure 5. Normalized autocorrelation functions for C(r = 12) and at /3 = 5.6 for the ensemble B 
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Figure 6. Normalized autocorrelation functions and integrated autocorrelation times for C(r = 12) (left) and 
Q2Q (right) at j3 = 5.8 for the ensemble D\ . 



relation time increases with the increasing size of the Wilson loop and also with the increasing 
smearing level. In the context of Wilson loop and Polyakov loop, SESAM collaboration has ob- 
served that geometrically extended observables suffer more from autocorrelation [19] with HMC 
algorithm. 

For hadronic observables the autocorrelations are quite small and since the number of mea- 
surements are not large the errors calculated from Eqs. (2.7) and (2.8) are quite large and mask the 
trends of the central values of autocorrelations. Our emphasis is on different trends of autocorrela- 
tions. To detect some trend of the central values of the autocorrelations for the hadronic observables 
we use a rough estimate of errors by single omission jackknife technique. In Fig. 1 1 integrated au- 
tocorrelation times for PP propagator with wall source for three different gaps (24,48,72) between 
measurements for the ensembles Bn a are presented. The data clearly exhibit the scaling properties 
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Figure 7. Normalized autocorrelation functions (left) and integrated autocorrelation times (right) for Pq at 
j8 = 5.6 for the ensembles Bn, and B$i,. 




Figure 8. Normalized autocorrelation functions (left) and integrated autocorrelation times (right) for 
^2o(4,4) at j3 = 5.6 for the ensembles B^ and B^i,- 

with the gaps. In Table 2 integrated autocorrelation times for pion (PP) and nucleon propagators 
with wall sources at a given time slice are presented. Clearly the integrated autocorrelation time 
decreases with increasing K (i.e. decreasing quark mass) both for pion and nucleon propagators. 
Similar observation was made by ALPHA collaboration in the case of Clover fermion [20]. The 
autocorrelation times of pion and nucleon propagators with point source and sink (not presented 
here) are smaller than the gap with which configurations are saved. 

For the determination of pion decay constants and PCAC quark mass, pion propagators other 
than PP are also needed. In Fig. 12 the integrated autocorrelation times for PP, AP, PA and AA 
correlators with wall source for the ensemble B?, a are presented. The propagators with A in the 
source are less correlated than P in the source. 

6 Improved estimation of T,„ f 

We have seen that the autocorrelations of different observables behave differently with the change 
in lattice spacing. As pointed out in [3], this behaviour is controlled by the coupling of different 
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Figure 9. Normalized autocorrelation functions and integrated autocorrelation times for Wilson loops with 
R = 1, T = 1 (left) and R = 4, T = 5 (right) at j8 = 5.8 for the ensemble D] . 
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Figure 10. Normalized autocorrelation functions and integrated autocorrelation times for Wilson loops with 
smearing levels = 5 (left) and smearing levels = 40 (right) at j3 =5.6 for the ensemble C%. 
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Table 2. Integrated autocorrelation times for pion (PP) and nucleon propagators with wall sources at j3 =5.6. 

observables with the slow modes of the transition matrix associated with Monte Carlo Markov 
chain. In this reference authors have proposed a method to quantify this coupling and estimate T, n( 
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Figure 11. Integrated autocorrelation times for PP propagator with wall source for three different gaps 
(24,48,72) between measurements for the ensemble Z?4„ at /3 = 5.6. 
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Figure 12. Integrated autocorrelation times for PP, AP, PA and AA correctors with wall source for the 
ensemble Z?3„ at j3 = 5.6. Measurements are done with a gap of 24 trajectories. 



more reliably. Following Ref. [3], an improved estimation of Ti nt can be determined as follows. 
Let T* be the best estimate of the dominant time constant. If for an observable G all relevant time 
scales are smaller or of the same order of T* then the upper bound of Tj nt 

<t = ~ +^ x T e {t)+A e {W u )%* (6.1) 
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Figure 13. Normalized autocorrelation function and effective autocorrelation time for Pq (left) £>2o (right) 
for the ensemble at j8 = 5.6. 

where Ag = max(r^(VK„),25r^(VK M )). W u is chosen where the autocorrelation is still significant. 
One possible estimation of T* is by measuring effective autocorrelation time, which is introduced 
in Ref. [3] as described below. Define effective exponential autocorrelation time 



exp 
V eff 



2 In 



r^(f/2) 



(6.2) 



z^ff which can be an estimate of T* is defined as, 



T ex P . 
l eff 



■ Maxff 



21n £^/2) 



(6.3) 



The estimation of X^FAG} requires good signal to noise ratio in the asymptotic region in a case by 
case basis which in turn requires very long Markov chain and is beyond the scope of the present 
work. 

However it is interesting to look at z^L{&) where reliable data is available and we present 
such an example in Fig. 13 (the jackknife technique is used to calculate the error of Z^L{&)). In 
Fig. 13 it appears that Q^ is coupling dominantly with slow mode, whereas Pq is coupling with 
more than one modes; nevertheless the slowest mode appearing in Pq is approximately the same as 
in Q% . This is reflected in the behaviour of T^y(^), which shows a single plateau for Q\ } , but for 
Pq, there is more than one plateau and the data is more noisy. Similar behaviour is observed in pure 
gauge theory in Ref. [3]. 

In improved estimation given in Eq. (6.1) central value of T,-, 1f gets modified. To check if 
this modification preserves the trend of autocorrelation of Q^o Wliri respect to quark mass, in Fig. 
14 we present the integrated autocorrelation times and their upper bounds {xf nt ) with rough errors 
estimated by jackknife method for topological susceptibilities (<2| ) at P = 5-6 (left) and at j3 = 5.8 
(right). At both lattice spacings, we find that both T, n ,(<22 ) and t" nt (Ql ) decrease as quark mass 
decreases. 
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Figure 14. Integrated autocorrelation times and their upper bounds (j" nt ) for topological susceptibilities 
(Qj ) at J3 = 5.6 (left) and at f5 = 5.8 (right). 

In conclusion, an extensive study of autocorrelation of several observables in lattice QCD 
with two degenerate flavours of naive Wilson fermion has shown that (1) at a given lattice spacing, 
autocorrelations of topological susceptibility and pion and nucleon propagators with wall source 
decrease with decreasing quark mass and autocorrelations of plaquette and Wilson loop do not in- 
crease with decreasing quark mass, (2) autocorrelation of topological susceptibility substantially 
increases with decreasing lattice spacing but autocorrelation of topological charge density corre- 
lator shows only mild increase and (3) increasing the size and the smearing level increase the 
autocorrelation of Wilson loop. 
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